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Abstract 

Adsorption of dimers is modelled using Random Sequential Adsorption algorithm. The inter- 
action between molecules is given by screened electrostatic potential. The paper focuses on the 
properties of adsorbed monolayers as well as the dependence of adsorption kinetics on interaction 
range. We designate random maximal coverage ratios, density autocorrelations and orientational 
ordering inside layers. Moreover the detailed analysis of adsorption kinetics are presented including 
discussion of Feder's law validity and new numerical method for modelling diffusion driven adsorp- 
tion. Results of numerical simulations are compared with experimental data obtained previously 
for insulin dimers. 
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I. INTRODUCTION 

Since its introduction by Feder [l| , Random Sequential Adsorption (RSA) became a well 
established method used for modelling of adsorption properties. Although at the beginning 
it was used mainly to model adsorption of simple spherical molecules, recent results shows 
that it could be effective also for quite complex structures like proteins {2, 3|. 

Most of the research effort focuses on adsorption of hard objects where geometry is the 
only factor affecting properties of obtained monolayers. On the other hand, adsorption is 
often induced by electrostatic interaction between adsorbate and collector, e.g. [4|. In such 
cases, the hard body interaction can still be sufficiently good approximation because the 
electrostatic forces are screened in a solution which makes them negligible. However in the 
general case, they should be taken into account to at least to find out the level of systematic 
error provided by such interactionless approximation. 

The purpose of presented paper is to extend previous investigations of dimers adsorption 
[st to include the case on non-negligible electrostatic repulsion. The paper focuses on fun- 
damental properties of dimer monolayers, such us maximal random coverage ratio, density 
autocorrelation and orientational ordering, as well as on adsorption kinetics. The additional 
aim is to develop robust numerical procedure to convert data obtained from RSA to values 
measured during typical adsorption experiment. 



II. MODEL 

A single dimmer molecule is assumed to consist of two identical, charged, spherical par- 
ticles. Electrostatic interaction between molecules originates from electrostatic potentials of 
those two particles. Typically, adsorbate molecules are dissolved in a polar solvent such us 
water, for example. In such conditions, the best approximation of electrostatic interaction 
is given by a screened electrostatic potential: 



2„2 



Uei{r) = exp 

4e^r 



r — a 



(1) 



where a and q stands for particle diameter and charge, respectively, with r denoting a 
distance between individual spheres belonging to different molecules. Qq is a constant pa- 
rameter and Sr characterises dielectric properties of a solvent. Range of the interaction is 



controlled by parameter commonly called as the Debye screening length. It is a key 
parameter for electrostatic interaction in electrolytes and can be calculated as: 



Le 



lekBT 
2eH 
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(2) 



where ks is the Boltzmann constant; T denotes temperature; e, elementary charge; and /, 
the ionic strength of electrolyte solution. Table [T] contains typical values of for the most 
common solutions. 



Concentration 

[M] 


1:1 

(KCl) 


1:2 [NaiSOi) 
2:1 (CaCh) 


2:2 

(NiSOi) 


1:3 (NasPOi) 
3:1 {AlCh) 


10-1 


0.9639 


0.5565 


0.4820 


0.3935 


10-2 


3.048 


1.759 


1.524 


1.244 


10-3 


9.639 


5.565 


4.819 


3.935 


10-4 


30.48 


17.59 


15.24 


12.44 


10-5 


96.39 


55.65 


48.19 


39.35 


10-6 


304.8 


175.9 


152.40 


124.4 



Table I. Le[nm] for typical electrolytes. Values were taken from 



Due to electrostatic repulsion, adsorption probability is not uniform over the whole col- 
lector. It depends on specific position and orientation of previously adsorbed, neighbouring 
molecules. Assuming thermal equilibrium of molecules in a bulk, this probability is similar 
to 

P{r) exp{-U{r)/kBT) (3) 

where f/(r) is a sum of Uel{r) of previously adsorbed particles. Figure [U shows adsorption 
probability ([3]) around a dimer for different Lg. 



A. Simulation details 



Adsorbed monolayers were generated using modified RSA algorithm. The procedure 
consists of the following steps: 

I: a new virtual dimer is randomly created. Its centre is set on the collector according to 
a uniform probability distribution and its orientation (the angle between X-axis and 
dimer axis) is uniformly chosen from [0,27r); 
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Figure 1. Probability of adsorption of a point-size elementary charge around a dimer molecule 
for different Debye screening length Lg = 0.1a, 1.0a and 5.0a. Parameters used: T = 293.15[i^], 



takes the role of distance unit. 

II. a: the virtual molecule undergoes overlapping test with its nearest neighbours; 

II. b: the total electrostatic potential U between the virtual molecule and previously ad- 
sorbed dimers is calculated using Eq.(IT]); 

III: if there is no overlap the virtual particle is added to the existing layer with a probability 
([3]). Otherwise the virtual dimer is removed and abandoned. 

The whole procedure is repeated for a specified number of times expressed using dimen- 
sionless time: 



where n is a number of algorithm iterations, Sm = is a coverage of a single molecule, 

and Sc is a collector's surface. The fundamental characteristic of an obtained layer is its 
coverage ratio defined as follows 



where Nm is a number of adsorbed molecules. 

The adsorption process simulation was performed for a squared collector of 200a-side 
size and was stopped at t = 10^. We did not used periodic boundary conditions, as it 
had been proved earlier, it does not have a significant influence on obtained layers [5;]. For 
each set of parameters, 20 to 100 independent simulations were performed. Parameters of 
electrostatic potential ([1]) were chosen to describe typical experimental conditions of water 



78 (as for water in a room temperature), and a = 9.29[nm\ which 




(4) 




(5) 



solutions. Therefore, the thermal energy ksT at a room temperature acts as energy unit 
and er = 78. The a = 9.29[nm] provides length scale typical to mid sized bio-molecules. 
Value of coefficient Qq = 6.78 [kBTa/e^] is fully determined by the above assumptions. 

III. RESULTS AND DISCUSSION 

A. Fundamental properties of adsorption monolayer 

Example coverages obtained in numerical simulations are shown in Figure Their main 




Figure 2. Typical adsorbed layers of dimers for three different electrostatic interaction ranges: 
Le = 0.1a, 1.0a and 5.0a. 

properties, such as maximal random coverage ratio, autocorrelation function and orienta- 
tional ordering, are analysed in the following sections. 



1. Adsorption ratio 

The easiest estimation of maximal random coverage ratio can be done by simple counting 
the number of adsorbed dimers. Figure |3] shows raw results taken directly from obtained 
data. 

However, such ratios are underestimated due to the finite time of a simulation; there is 
no guarantee that all free places on a collector have been filled. To deal with this, the model 
of RSA kinetics have to be used. The common choice here is the Feder's law 17- 91, which 

n 

is valid for a wide range of adsorbate molecule's shapes [10[ and also proved to be valid for 
the case of hard-core dimers RSA [sl: 

9ma. - e{t) ~ r^/^ (6) 
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Figure 3. Adsorption ratio dependence on electrostatic interaction range Lg. The ratio is a mean 
value of 20 to 100 simulations using a collector of 200a side size. The time of a single simulation 
was 10^. 

where t is dimensionless time @ and d is a dimension of a collector. Here, d = 2. Plots 
illustrating relation (Q are presented in Figure |H They confirm that the numerical data 
obey the Feder's law. Moreover, at the limit of t — )■ oo (t^^/^ — > 0), the maximal random 
coverage ratios are higher by about 1-2% than the values obtained directly from data as in 
Figure [3] (see Table HTl) . 
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Figure 4. Coverage ratio 9 dependence on t~^^'^ for different electrostatic interaction range. Lg G 
{0.1, 0.2, 0.5, 1.0, 2.0, 5.0}. Black dots are taken from simulation and red line is a fit dS]). The 6max 
is reached when t^^^'^ = 0. 



It is clear that the maximal random coverage ratio defined using ([5]) and decreases 
with the growth of the Debye screening length due to the electrostatic repulsion. On 



the other hand, those results can also be interpreted as an increase in the effective molecule 
size with 6max = 0.547 being constant. The effective molecule size can be also estimated 
analytically using Eq.([n]): 



Seff{Le 



1 - exp[-Uei{Le,r) / ksT] if exp [~Uei{Le,r) /ksT] > a, 







otherwise. 



(7) 



Seff counts the area where adsorption probability is higher than a. For a — > 0+ the original 
Sm = 110^/2 is reproduced. Here, a = 0.02 was used to get the best fit for the numerical 
data. As it has been shown in TabJIT| the above analytic approximation matches the data 
within 10% error margin. 
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a = 0.02) [a^] 


0.0 


0.547 


0.002 




1.57 


1.57 


0.1 


0.51008 


0.00005 


0.9994 


1.69 


1.73 


0.2 


0.48535 


0.00014 


0.9975 


1.77 


1.87 


0.5 


0.43538 


0.00017 


0.9983 


1.97 


2.25 


1.0 


0.37185 


0.00020 


0.9987 


2.31 


2.81 


2.0 


0.27660 


0.00019 


0.9990 


3.11 


3.73 


5.0 


0.14623 


0.00015 


0.9968 


5.87 


5.44 



Table II. Maximal random coverage ratio and effective molecule area for different values of after 
applying the Feder's law correction. Correlation coefficient describes agreement between numerical 
data and relation ([6]). 1 or —1 means perfect fit. 



2. Autocorrelations 

Autocorrelation function G(r), also known as two point correlation function, is defined 
as a mean density of adsorbed dimers at a given distance from the centre of one adsorbed 
molecule. Here, the density is normalised to the overall mean density in a covering layer. 
Therefore, the autocorrelation function presented in Figure [5] approaches 1.0 with growing 
distance. Autocorrelation behaviour is typical. For small electrostatic interaction range, the 
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Figure 5. Spatial autocorrelation function. It was obtained for a collector size 200a x 200a with 
different interaction range Lg. All functions are normalised to be equal to 1 at infinity. 

rapid grow starts at r = 2. As expected, an increase in Lg makes this grow slower. For 
Le = 5.0, there are almost no neighbours in the distance closer than r = 3. The maximum of 
G{r) corresponds to the most probable distance between molecules. With the grow of Lg it 
moves to the right and becomes wider, but only for high values. The data are presented 
in Table IIIII For < 1.0, the maximum is followed by the minimum, which corresponds 





characteristic distance 


dispersion 


0.1 


3.68 


1.215 


0.2 


3.80 


1.134 


0.5 


4.00 


1.215 


1.0 


4.29 


1.296 


2.0 


4.98 


1.215 


5.0 


6.84 


1.539 



Table III. Characteristic distance between neighbouring dimers and it's dispersion (FWHM) de- 
pendence on interaction range Lg. 



1.5 



O 

0.5 

n, 



to the excluded volume effect around nearest neighbours. Careful reader may also observe 
here second smaller maximum. 

The phenomenological approach shows that autocorrelation function for spheres is char- 
acterised by a log singularity at a kiss limit (r — > 2+) and superexponential decay for r — t- cxo 
This results were also confirmed for non-interacting dimers |5|. However, screened 



electrostatic potential makes it difficult to precisely determine the kiss limit. Moreover, the 
long range limit is also affected because the electrostatic interaction, obviously, vanishes not 
as fast as hard-core interaction. Therefore in this case, the superexponential decay cannot 
be expected. 

3. Order parameter 

Dimers have a non-uniform shape and therefore some orientational ordering may occur 
in adsorbed layers. This phenomenon is important from a practical point of view although 
it affects important mechanical, electrochemical and optical properties of a formed layer. 
To determine whether any orientational order appears in a coverage, measurement of such 
an ordering is needed. It can be based on the function S{a), which is the sum of squared 
scalar products between molecules orientation and a given angle: 

1 ^ 

S{a) = —'^{xi ■ cosa + Hi ■ sina) , (8) 

^ 1=1 

where {xi,yi) denotes a unit vector along the direction of i-th molecule in a layer. In the 
case of all molecules oriented along the specified direction aQ, S{a) will have a maximum 
S{ao) = S{ao + vr) = 1 and a minimum S{aQ + 7r/2) = S{ao + 37r/2) = 0. When there is 
no order, there will be no extremes and S{a) ~ 0.5. The mean orientation is given by the 
maximum of S{a) function and can be estimated from: 

tan2ao=wv '^^_Viv > (9) 

and the order parameter can be defined as S* = 5'(q;o). It is worth to notice that the above 
equation is satisfied by both ao, and ao + t{/2; so it should be determined which one is 
a maximum and which is a minimum. Dependence between 5* on interaction range L^. is 
shown in Figure [61 All the values are close to S" = 0.5, which means the orientations of 
particles are randomly distributed and there is no significant global ordering. However, 
there is small but noticeable growth of 5* with the growth of interaction range. We have 
checked that this effect is more evident for smaller collectors sizes. It is probably induced 
by adsorption conditions near collector boundaries. With the grow of interaction range the 
density of adsorbed molecules decreases and the influence of boundaries propagates on larger 
distances, which explains the growth of global ordering. 
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0.516 




Figure 6. Order parameter for different interaction range. The presented value is a mean taken 
from 20 to 100 simulations depending on collector size. The time of single simulation was 10^. 



To complete this picture, we measured also two-point correlation functions of the orien- 
tational order. It allows the measurement of local orientational ordering. Results in Fig. 
[7] present a mean value of the scalar product of main particle axes separated by a given 
distance. When interaction range is small (Le < 0.5) the nearest neighbours slightly prefer 




Figure 7. Order autocorrelation function for different Debye screening lengths L^. 

parallel alignment; whereas when Le > 1.0, the S{r) drops below 0.5 which means slight 
tendency to perpendicular alignment. For larger distances, all plots quite quickly approach 
0.5, which is accordance with previous observations and means lack of global ordering. 
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B. Adsorption kinetics 



The growth of covering layer for homogeneous systems can be described by one dimen- 
sional differential equation: 

j^e{t) = k-AFS{e)-c{Q,t), (10) 

where ASF {9), known also as Available Surface Function, is a ratio of uncovered adsorption- 
ready surface to the whole collector area; k is an adsorption reaction constant; and c(0,t) 
is a molecules concentration near the collector surface. The RSA kinetics assumes that the 
only factor affecting the speed of adsorption layer growth is decreasing probability of finding 
large enough, unallocated place on the collector surface. From a physical point of view, it 
is equivalent to the assumption that concentration c(0, t) of molecules in a solution near 
the layer is constant. However, the value of c(0,t) is generally a result of different physical 
processes which move particles from a bulk to a surface neighbourhood. Their details depend 
on an experiment preparation and its environment. On the other hand, the most common 
process of transport in the majority of experiments is diffusion: 

dc{x,t) _ d'^c{x,t) 

dt dx^ ' ^ ' 

where D is a diffusion constant. The boundary conditions for the above equation are defined 
by constant bulk concentration and continuity of particles flux at the surface: 

lim cix, t) = Coo (12) 

- = k ■ ASF{9) ■ c(0, t) (13) 

ox 

Equations ([T^ and f lT^ are known as mixed boundary conditions or Robin boundary con- 
ditions. Due to nonlinearity of the ASF (9), the equation set f|T0l) - f|T3l) can be solved only 
numerically. The details of numerical procedure used here are described in Appendix lAl 



1. Measurement of the ASF{0) 

The first step to determine real time adsorption kinetics 9{t) governed by Eqs. fllip - fll3p 
is to find out Available Surface Function ASF{9). It can be done by estimating ratio of 
successful to all RSA attempts for a given coverage. The total number of attempts was set 
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typically as 10^; however, if all of the attempts failed the simulation was performed till the 
first successful attempt. Obtained results are presented in Figure [HI The ASF{9) behaviour 




Figure 8. Available Surface Function for different range of electrostatic interaction. Inset shows 
the same data in a log scale. Solid lines correspond to fit given by Eqll4[ 

is typical. In the small coverages limit ASF{6) can be approximated by a linear or square 
function of 6. For jamming limit, the exponential decrease is observed. To interpolate 
ASF{6) in the whole range of 6, we have chosen a fit used previously for ellipsoids and 
sphero cylinders 



where 6 



ASF{e) = (1 + aiO + 02^^ + a3^^)(l - Of, (14) 
The values of fitted parameters are presented in Table HVl Relation ([H]) 





ai 


02 


as 


correlation coefficient 


0.1 


1.19711 


1.97556 


-4.20903 


0.999993 


0.2 


1.06048 


1.53124 


-5.09405 


0.999996 


0.5 


0.69920 


-1.50719 


-2.64705 


0.999989 


1.0 


-0.01671 


-5.37355 


4.41001 


0.999976 


2.0 


-1.67599 


-4.50836 


8.04173 


0.999977 


5.0 


-4.40441 


4.43515 


0.983194 


0.999841 



Table IV. Fitted parameters according to EqlTll 

gives the best fit when interaction range is small. For large Lg the ASF quickly approaches 
exponential function rather than polynomial one, which can be clearly seen when using 
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logarithmic scale. In general, the fit given by (fT^ breaks when ASF becomes exponential. 
For Le = 0.1, it takes place at ^ ~ 0.7, whereas for Lg = 5.0 at ^ ^ 0.3. 



2. Real time kinetics 

Numerical results of dimers adsorption were compared with experimental data of insulin 

I i 

adsorption obtained by MoUmann et al. [12]. Therefore, we set the diffusion coefficient 
D = 100 [/im^/s] and dimer surface Sm = 8.75 ■ 10~^ [/i^]. The bulk concentration was 
10^'^ [mg/ 171"^]. The value of adsorption reaction rate had not been specified explicitly; there- 
fore, we used k = 1 [jj,/ s] in our calculation. Results are presented in Figure[9l The units used 
here are commonly used by experimenters. Presented plots drawn for different electrostatic 




1000 2000 

time [s] 



3000 



Figure 9. Real time kinetics of dimers adsorption. Different lines corresponds to different Lg values. 
Dots are taken from experiment jl^ . 



interaction range are similar. The most noticeable difference is connected with maximum 
coverage ratio, which strongly depends on interaction range (see fFigj3]and Tab HTj) . How- 
ever, it can be noticed that the saturation occurs faster with the growth of interaction range, 
which is a direct consequence of the ASF{6) shape. The difference between simulation and 
experimental data result mainly from different values of the maximal random coverage ra- 
tios. The causes of this effect were already discussed earlier jst . Apart from this, the shapes 
of numerically and experimentally determined kinetics are similar. 

13 



IV. SUMMARY 



The random maximal coverages ratio for electrostatically interacting dimers was measured 
using the RSA method. Repulsive interactions lowers the maximum possible coverages ratio 
the more the larger interaction range Le is, eg. 9max = 0.510 for Le = 0.1 and 9max = 0.146 
for Le = 5.0; however, it does not affect the validity of the Feder's law. The measurement of 
autocorrelation function allows to find out the characteristic distance between neighbouring 
molecules and its value changes from 3.68 for = 0.1 to 6.84 for = 5.0. The spontaneous 
orientational ordering is imperceptible at the global scale but it can occur locally and can 
be exaggerated by the presence of boundaries. The adsorption kinetics depends mainly on 
maximal coverage ratio but it can be observed that saturation value is approached faster for 
larger Lg. 

This work was supported by MNiSW/0013/H03/2010/70. 



Appendix A: Evaluating the kinetics of adsorption 



'here are known some algorithms for calculating adsorption kinetics in case of adsorption 
1^ : however, we want to introduce here another scheme based only on diffusion equations 
(fT0|) - (fT3!l . which can be easily applied to any chemical reaction. Therefore, the diffusion 
equation is solved using standard Crank-Nicholson method; however, the introduction of 
boundary conditions requires additional operations after each step of the algorithm, which 
are described below: 

i) : the concentration having the largest x coordinate is set to Coo- 

ii) : the adsorption probability is given by: 

p{t) = k ■ ASF{e) ■ At/ Ax. (Al) 
It has been proved that the above relation fully satisfies Robin boundary conditions 



15| . Following relations are direct consequences of the above. The number of 



particles near the surface is equal to 



n{t) = c{0,t) ■ Se- Ax. (A2) 
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Thus, the number of adsorbed molecules is: 



n^(t) = n{t) ■ p{t) = c(0, t)-S,-k- ASF{e) ■ At. (A3) 

Ax and At are space and time discretization steps used in Crank-Nicholson algorithm. 
The coverage increase corresponding to the above number ha is given by: 

Ae{t) = nA{t) ■ Sm/Sc = c(0, t) ■ kASF{9) ■ At ■ (A4) 

so the concentration near surface changes according to 

c(0, t + At) = c(0, t) - nA{t)/{Sc ■ Ax) = c(0, t)[l - k ■ ASF{e) ■ At/ Ax] (A5) 

To use the above algorithm, the diffusion constant D, Available Surface Function ASF{9) 
and reaction rate k is needed. The diffusion constant depends on adsorbed molecule prop- 
erties (shape, mass) as well as solvent properties and it can be measured experimentally. 
ASF{9) can be easily determined from RSA simulation as a ratio of successful to all adsorp- 
tion attempts on a given layer. The adsorption reaction rate k can be expressed in terms 
of physical parameters characterising the system such as the particle diffusion coefficient, 
specific energy distribution, depth equilibrium state, and height of the adsorption energy 
barrier jsl, Q. For barrier- less adsorption regime, the adsorption constant is given by 

k = — — (A6) 
where 6a is the thickens of adsorption boundary and 6m is the minimum distance between 



molecule in bulk phase and surface [16!]. Reaction rate can also be measured experimentally 
or used as a parameter which can be fitted to obtain the best possible agreement between 
experiment and numerical simulation. 
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